# THIS CODES PRODUCES THE DATA FOR TABLE X IN THE PAPER
# "Where's the Bias", by Cesar Zucco Jr.
#
# It simulates the seat share obtained by one party against a unified opposition with varrying magnitude
# More details on the simulation write up. 
#
# The expected vote share of the party is set at m=43 (from Navia)
# The "for var." stipulates the variance levels that will be simulated (currently 3 levels)
# The function transf computes the appropirate alpha and beta of the beta distribution given m and the variance level
# Once the beta distribution for that variance level is defined, 
#        (1) I allow for different magintudes in the "for mag"" loop
#        (2) I stipulate the apporiate number and boundary of brackets to obtain 0,1,...,mag seats 
#        (3) I compute the c.d.f beta for each bracket in the "for i" loop
#
# Bottomline: Magnitude = 2 was not the best system a priori



#SIMULATION
par(mfrow=c(3,1),mar=c(4,4,1,1),ask=FALSE) #mar(bottom, left, top, right)
m = 0.43
marker <- 1
vars <- c(0.0025,0.01,0.0225)
overall <- matrix(nrow=10,ncol=6)
colnames(overall) <- c("",sqrt(vars[1]),"",sqrt(vars[2]),"",sqrt(vars[3]))
for(var. in vars){     #Determines that sd levels 0.05, 0.1 and 0.15 respectively
                                      
    transf <- function(a,x,v) {         #this function was obtained from the first and second moments of the binomial - see explanation
        (a*(a+1))/((a+((a*(1-x))/x))*(a+((a*(1-x))/x)+1)) - x^2 -v          
    }

    temp <- uniroot(transf
        , interval=c(0.01,100000)
        , tol = 0.0001
        , x=m           #x is the mean vote share across districts, that varies in each "for m" loop
        , v=var.)       #v is the standard deviation, that varies in each of the levels set at the "for var." loop                       
    new.a <- temp$root
    new.b <- (new.a*(1-m))/m


expec.seat.share <- matrix(nrow=10,ncol=2)
colnames(expec.seat.share) <- c("Mag","S.Share")
for(mag in c(1:10)){
brackets <- 1+mag
each.bracket <- matrix(nrow=mag+2,ncol=7)
colnames(each.bracket) <- c("mag","bracket","cut","cum.pr","prob","seats","seat.share")
bracket <- 0
bounds <- seq(0,1,1/(1+mag))
for(i in seq(0,1,1/(1+mag))){
    if (i == 0)  { 
       cum.bracket <- 0
       seats <- 0
       pr.bracket <- 0
       s.share <- 0
   } else {    
       cum.bracket <- pbeta(i,new.a,new.b)
       seats <- bracket - 1
       s.share <- seats/mag
       pr.bracket <- cum.bracket - each.bracket[bracket,4] 
   }    
   each.bracket[bracket+1,] <- c(mag,bracket,i,cum.bracket,pr.bracket,seats,s.share) 
   bracket <- bracket + 1   
   }
expec.seat.share[mag,] <- c(mag,sum(each.bracket[,5]*each.bracket[,7]))
}
overall[,c(marker,marker+1)] <- expec.seat.share
marker = marker + 2
}
cat("\n","\n","Table X in the paper","\n","\n")
print(round(overall[ ,c(1,2,4,6)],3))
library(xtable)
cat("\n","\n","Table for LATEX","\n","\n")
print(xtable(round(100*overall[,c(2,4,6)],2)))
